#!/usr/bin/python
# -*- coding: utf-8 -*-

# Copyright (c) 2011
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
#
# Author: Jesus Carrero <j.o.carrero@gmail.com>
# Mountain View, CA

from scipy.special import gammaln
from Polynomials import _jacobiP
import numpy as np
import fileinput


class Quadrature(object):
    """ Base class for all quadrature formulas. """

    __slots__ = ['m_weights', 'm_points', 'm_order', 'm_data']

    def __init__(self, points, weights, order):
        """ Base class to handle all quadrature formulas. """

        self.m_weights = weights
        self.m_points = points
        self.m_order = order
        self.m_data = zip(points, weights)

    def __call__(self, func):
        """Integrate the callable f with respect to the given quadrature rule.
        """

        return sum(w * func(x) for (x, w) in self.m_data)

    def get_order(self):
        """ return quadrature order. """

        return self.m_order

    def get_weights(self):
        """ quad weights. """

        return self.m_weights

    def get_qpoints(self):
        """ quadrature points. """

        return self.m_points


class Simplex2DQuadrature(Quadrature):
    """ 2D quadrature formulas for triangles. """

    def __init__(self, order):
        """ Load quadrature from tables. """

        (nodes, weights) = self.load_quads(order)
        Quadrature.__init__(self, nodes, weights, order)

    @staticmethod
    def load_quads(order):
        """ load a quadrature formula of the given order from a table """

        if 1 > order | order > 7:
            print 'Choose an quad order between 1 and 7'
            return

        tables = fileinput.input('Tables.txt')
        for line in tables:
            row = line.split()
            if len(row) == 0:
                continue
            if line[0] == '#':
                continue

            l = [int(x) for x in row]
            if l[1] == order:
                points = np.zeros((l[2], 2), 'float64')
                weights = np.zeros((l[2], 1), 'float64')
                for s in range(0, l[2]):
                    row = tables.readline().split()
                    points[s, :] = [float(x) for x in (row[0], row[1])]
                    weights[s] = float(row[2])

                break
            else:
                for s in range(0, l[2]):
                    row = tables.readline().split()

        fileinput.close()
        return (points, weights)


class JacobiGaussQuadrature(Quadrature):

    """An M{N}th order Gauss quadrature associated with the Jacobi
    polynomials of type M{(alpha,beta) > -1}

    C{alpha} and C{beta} may not be -0.5.

    Integrates on the interval (-1,1).
    """

    def __init__(self, alpha, beta, order):
        """ Look at the book: Spetral Methods. """

        (points, weights) = self.compute_weights_and_nodes(order, alpha,
                beta)
        Quadrature.__init__(self, points, weights, 2 * order - 1)

    @staticmethod
    def compute_weights_and_nodes(order, alpha, beta):
        """Return (nodes, weights) for an n-th order Gauss quadrature
        with the Jacobi polynomials of type (alpha, beta).
        """

        nodes = np.zeros((order + 1, 1), 'float')
        weights = np.zeros((order + 1, 1), 'float')

        for k in np.linspace(0, order, order + 1, endpoint=True):
            r = -np.cos(np.pi * (2.0 * k + 1.0E0) / (2.0 * order))
            if 0 < k:
                r = (r + nodes[k - 1]) / 2.0
            epsilon = 1.0e-16
            delta = 1.0E0
            while epsilon < abs(delta):
                s = np.sum(np.array([1 / (r - nodes[m]) for m in
                           np.linspace(0, k, k, False)]))
                (y, dy, d2y) = _jacobiP(order + 1, alpha, beta, r)
                delta = -y / (dy - y * s)
                r += delta
            nodes[k] = r

        ab1 = alpha + beta + 1.0E0
        ao1 = alpha + order + 1.0E0
        bo1 = beta + order + 1.0E0
        abo1 = ab1 + order
        aux = ab1 * np.log(2.0) + gammaln(ao1) + gammaln(bo1) \
            - gammaln(abo1)
        aux -= np.sum(np.array([np.log(i) for i in np.linspace(2,
                      order, order - 1, True)], 'float'))
        aux = np.exp(aux)

        for i in np.linspace(0, order, order + 1, True):
            (y, dy, d2y) = _jacobiP(order + 1, alpha, beta, nodes[i])
            weights[i] = aux / ((1.0E0 - nodes[i] ** 2) * dy ** 2)

        return (nodes, weights)


class LegendreGaussQuadrature(JacobiGaussQuadrature):

    """M{N}th order Gauss quadrature associated with the Legendre polynomials.
    """

    def __init__(self, order):
        JacobiGaussQuadrature.__init__(self, 0, 0, 2 * order - 1)


class TransformedQuadrature1D(Quadrature):

    """A quadrature rule on an arbitrary interval M{(a,b)}. """

    def __init__(self, quad, left, right):
        """Transform a given quadrature rule `quad' onto an arbitrary
        interval (left, right).
        """

        self.m_left = left
        self.m_right = right

        length = right - left
        assert length > 0
        half_length = length / 2
        Quadrature.__init__(self, [left + (p + 1) / 2 * length for p in
                            quad.m_points], [w * half_length for w in
                            quad.m_weights])


if __name__ == '__main__':
    quad = LegendreGaussQuadrature(4)

    print ' ( INFO ) Testing 1D quadrature points '
    print ' ( INFO ) This is the quad formula for n = 10'
    print """ [ qpoints   |   weights ] """
    print np.c_[quad.get_qpoints(), quad.get_weights()]
    print """ scaled to [0, 1] """
    print np.c_[(1.0E0 + quad.m_points) / 2.0, quad.m_weights]

    # quad2 = Simplex2DQuadrature( 2 )
